## [1] "Run Completed at 2016-07-29 14:53:50"
#Load in data from Observed2m_Generate.Rmd
load("Observed.RData")
#load("ObservedModel.RData")
newModel<-T
#source functions
source("Bayesian/BayesFunctions.R")

View Raw Data

ggplot(indat,aes(x=Traitmatch,y=Yobs,col=as.factor(R))) + facet_wrap(~Hummingbird,ncol=4,scales="free") + geom_point() + geom_smooth(method = "glm",method.args=list(family="poisson")) + scale_color_manual("Resource Availability",values=c("black",'blue','red'))

ggplot(indat[indat$Yobs==1,],aes(x=Traitmatch,fill=as.factor(R))) + facet_wrap(~Hummingbird,ncol=4,scales="free") + geom_density(alpha=.7) + scale_fill_manual("Resource Availability",values=c("black",'blue','red'))

m<-glmer(data=indat[,],Yobs~Traitmatch*R+(1|Hummingbird),family="poisson")
summary(m)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Yobs ~ Traitmatch * R + (1 | Hummingbird)
##    Data: indat[, ]
## 
##      AIC      BIC   logLik deviance df.resid 
##   6372.1   6414.8  -3179.0   6358.1     3314 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0292 -0.6257 -0.3757 -0.2129 26.7497 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  Hummingbird (Intercept) 1.536    1.24    
## Number of obs: 3321, groups:  Hummingbird, 19
## 
## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -0.610860   0.296018  -2.064   0.0391 *  
## Traitmatch         -0.032380   0.004328  -7.482 7.33e-14 ***
## RMedium            -0.171362   0.094937  -1.805   0.0711 .  
## RHigh              -0.206707   0.088165  -2.345   0.0191 *  
## Traitmatch:RMedium -0.003301   0.006738  -0.490   0.6242    
## Traitmatch:RHigh   -0.002527   0.006241  -0.405   0.6856    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmtc RMedim RHigh  Trt:RM
## Traitmatch  -0.162                            
## RMedium     -0.152  0.484                     
## RHigh       -0.165  0.532  0.507              
## Trtmtch:RMd  0.108 -0.646 -0.735 -0.347       
## Trtmtch:RHg  0.119 -0.706 -0.347 -0.737  0.461
## convergence code: 0
## Model failed to converge with max|grad| = 0.00161504 (tol = 0.001, component 1)
indat$pred<-predict(m,type="response")
ggplot(indat,aes(x=Traitmatch)) + geom_point(aes(y=Yobs)) + geom_line(aes(y=pred,col=as.factor(R)))

Hierarchical Occupancy Model

For hummingbird species i feeding on plant species j observed at time k and sampling event observed by transect (YTransect) or camera (YCamera)

Observation Model:

\[ YCamera_{i,j,k,d} \sim Binomial(N_{i,j,k},\omega_{Camera}) \] \[ \omega_{Camera} <- \phi_{Camera} * EffortCamera_k \]

Process Model:

\[ N_{i,j,k} \sim Poisson(\lambda_{i,j,k}) \] \[ log(\lambda_{i,j,k}) = \alpha_i + \beta_{1,i} * Traitmatch_{i,j} + \beta_{2,i} *Resources_k + \beta_{3,i} * Traitmatch_{i,j} * Resources_k \]

Priors

\[ \phi_{Camera} \sim U(0,1) \] \[\alpha_i \sim Normal(\alpha_\mu,\alpha_\tau)\] \[\beta_{1,i} \sim Normal(\mu_{\beta_1,\tau_{beta_1}})\] \[\beta_{2,i} \sim Normal(\mu_{\beta_2,\tau_{beta_2}})\] \[\beta_{3,i} \sim Normal(\mu_{\beta_3,\tau_{beta_3}})\]

Group Level Means \[ \mu_\alpha \sim Normal(0,0.0001)\] \[\mu_{\beta_1} \sim Normal(0,0.0001)\] \[\mu_{\beta_2} \sim Normal(0,0.0001)\] \[\mu_{\beta_3} \sim Normal(0,0.0001)\]

Group Level Variance \[\tau_{\alpha} \sim Uniform(0,1000)\] \[\tau_{\beta_1} \sim Uniform(0,1000)\] \[\tau_{\beta_2} \sim Uniform(0,1000)\] \[\tau_{\beta_3} \sim Uniform(0,1000)\]

#Source model
source("Bayesian/NmixturePoissonRagged2m.R")

#print model
writeLines(readLines("Bayesian/NmixturePoissonRagged2m.R"))
## 
## sink("Bayesian/NmixturePoissonRagged2m.jags")
## 
## cat("
##     model {
##     #Compute true state for each pair of birds and plants
##     for (i in 1:Birds){
##     for (j in 1:Plants){
##     for (k in 1:Times){
##     
##     #Process Model
##     log(rho[i,j,k])<-alpha[i] + beta1[i] * Traitmatch[i,j] + beta2[i] * resources[i,j,k] + beta3[i] * resources[i,j,k] * Traitmatch[i,j]
##     
##     #True State
##     S[i,j,k] ~ dpois(rho[i,j,k])
##     }
##     }
##     }
##     
##     #Observation Model
##     for (x in 1:Nobs){
##     
##     #Observation Process for cameras
##     Yobs_camera[x] ~ dbin(dcam[Bird[x]],S[Bird[x],Plant[x],Time[x]])    
## 
##     #     #Assess Model Fit - Posterior predictive check
##     # 
##     #     #Fit discrepancy statistics
##          #Camera
##          eval_cam[x]<-dcam[Bird[x]]*S[Bird[x],Plant[x],Time[x]]
##          E[x]<-pow((Yobs_camera[x]-eval_cam[x]),2)/(eval_cam[x]+0.5)
##     #     
##          ynew_cam[x]~dbin(dcam[Bird[x]], S[Bird[x],Plant[x],Time[x]])
##          E.new[x]<-pow((ynew_cam[x]-eval_cam[x]),2)/(eval_cam[x]+0.5)
##     }
##     
##     #Priors
##     #Observation model
##     #Detect priors, logit transformed - Following lunn 2012 p85
##     
##     for(x in 1:Birds){
##     #For Cameras
##     logit(dcam[x])<-detect[x]
##     detect[x] ~ dnorm(dprior_cam,tau_dcam)
##     }
##     
##     #Detection group prior
##     dprior_cam ~ dnorm(0,0.386)
## 
##     #Group effect detect camera
##     tau_dcam ~  dt(0,1,1)I(0,)
##     sigma_dcam<-pow(1/tau_dcam,.5)
##     
##     #Process Model
##     #Species level priors
##     for (i in 1:Birds){
##     
##     #Intercept
##     #logit prior, then transform for plotting
##     alpha[i] ~ dnorm(alpha_mu,alpha_tau)
##     
##     #Traits slope 
##     beta1[i] ~ dnorm(beta1_mu,beta1_tau)    
##     
##     #Plant slope
##     beta2[i] ~ dnorm(beta2_mu,beta2_tau)    
##     
##     #Interaction slope
##     beta3[i] ~ dnorm(beta3_mu,beta3_tau)    
##     }
##     
##     #Group process priors
##     
##     #Intercept 
##     alpha_mu ~ dnorm(0,0.386)
##     alpha_tau ~ dt(0,1,1)I(0,)
##     alpha_sigma<-pow(1/alpha_tau,0.5) 
##     
##     #Trait
##     beta1_mu~dnorm(0,0.386)
##     beta1_tau ~ dt(0,1,1)I(0,)
##     beta1_sigma<-pow(1/beta1_tau,0.5)
##     
##     #Resources
##     beta2_mu~dnorm(0,0.386)
##     beta2_tau ~ dt(0,1,1)I(0,)
##     beta2_sigma<-pow(1/beta2_tau,0.5)
##     
##     #Interaction
##     beta3_mu~dnorm(0,0.386)
##     beta3_tau ~ dt(0,1,1)I(0,)
##     beta3_sigma<-pow(1/beta3_tau,0.5)
## 
##     #derived posterior check
##     fit<-sum(E[]) #Discrepancy for the observed data
##     fitnew<-sum(E.new[])
##     
##     }
##     ",fill=TRUE)
## 
## sink()
#Inits
InitStage <- function(){
  #A blank Y matrix - all present
  initY<-array(dim=c(Birds,Plants,Times),22)
  initB<-rep(0.5,Birds)
list(S=initY,detect=initB)}

#Parameters to track
ParsStage <- c("alpha","beta1","beta2","beta3","alpha_mu","alpha_sigma","beta1_sigma","beta1_mu","beta2_mu","beta2_sigma","beta3_mu","beta3_sigma","dcam","dprior_cam","E","fit","fitnew")


#Jags Data
Dat<-list(
  Yobs_camera = indat$Camera,
  Birds=max(indat$jBird),
  Bird=indat$jBird,
  Plant=indat$jPlant,
  Time=indat$jID,
  Plants=max(indat$jPlant),
  Times=max(indat$jID),
  resources=resourceMatrix,
  Nobs=nrow(indat),
  Traitmatch=jTraitmatch)

  #MCMC options
  if(newModel){
    system.time(
      m2<-jags.parallel(data=Dat,parameters.to.save =ParsStage,inits=InitStage,model.file="Bayesian/NmixturePoissonRagged2m.jags",n.thin=4,n.iter=200000,n.burnin=200000-2000,n.chains=2,DIC=F)
      )
  }
##     user   system  elapsed 
##   11.079    0.422 5267.675
#recompile if needed
runs<-100000

recompile(m2)

if(!newModel){
  system.time(m2<-update(m2,n.iter=runs,n.burnin=runs-2000,n.thin=2))  
}
#extract par to data.frame
pars_detect<-extract_par(m2,data=indat,Bird="jBird",Plant="jPlant")

Assess Convergence

###Chains
ggplot(pars_detect[pars_detect$par %in% c("alpha","beta1","beta2","beta3"),],aes(x=Draw,y=estimate,col=as.factor(Chain))) + geom_line() + facet_grid(par~species,scale="free") + theme_bw() + labs(col="Chain") + ggtitle("Species Level")

ggplot(pars_detect[pars_detect$par %in% c("dcam"),],aes(x=Draw,y=estimate,col=as.factor(Chain))) + geom_line() + facet_grid(species~par,scale="free") + theme_bw() + labs(col="Chain") + ggtitle("Species Level")

ggplot(pars_detect[pars_detect$par %in% c("dprior_cam","beta1_mu","beta1_sigma","beta2_mu","beta2_sigma","beta3_mu","beta3_sigma","alpha_mu","alpha_sigma"),],aes(x=Draw,y=estimate,col=as.factor(Chain))) + geom_line() + theme_bw() + labs(col="Chain") + ggtitle("Group Level Parameters") + facet_wrap(~par,scales="free")

Posteriors

###Posterior Distributions
ggplot(pars_detect[pars_detect$par %in% c("alpha","beta1","beta2","beta3"),],aes(x=estimate)) + geom_histogram(position='identity') +  facet_grid(species~par,scales="free") + theme_bw() + ggtitle("Species Parameters")

#Detection figure
pars_detect<-merge(pars_detect,jagsIndexBird,by.x="species",by.y="jBird",all=T)
ggplot(pars_detect[pars_detect$par %in% "dcam",],aes(x=par,y=estimate)) + geom_violin(fill='black') + theme_bw() + ggtitle("Detection Probability") + facet_wrap(~Hummingbird)

Detection Table

detecttable<-group_by(pars_detect,Hummingbird,par) %>% filter(par %in% c('dcam')) %>% summarize(mean=mean(estimate),lower=quantile(estimate,0.05),upper=quantile(estimate,0.95))
detecttable
## Source: local data frame [19 x 5]
## Groups: Hummingbird [?]
## 
##                  Hummingbird    par      mean      lower     upper
##                       (fctr) (fctr)     (dbl)      (dbl)     (dbl)
## 1             Andean Emerald   dcam 0.2756901 0.12795174 0.5023129
## 2         Booted Racket-tail   dcam 0.2552547 0.16037808 0.3802274
## 3                 Brown Inca   dcam 0.2240020 0.12590289 0.3608807
## 4        Buff-tailed Coronet   dcam 0.2284291 0.11936447 0.3819514
## 5              Collared Inca   dcam 0.2559656 0.13416553 0.4183749
## 6          Crowned Woodnymph   dcam 0.2782735 0.14273788 0.4410311
## 7    Fawn-breasted Brilliant   dcam 0.2703968 0.10042306 0.4995664
## 8          Gorgeted Sunangel   dcam 0.5285593 0.25151338 0.7797808
## 9    Green-crowned Brilliant   dcam 0.2306932 0.07970214 0.4315912
## 10   Green-fronted Lancebill   dcam 0.3600551 0.23995363 0.4987622
## 11             Hoary Puffleg   dcam 0.3124265 0.15226076 0.5507799
## 12    Purple-bibbed Whitetip   dcam 0.2735183 0.12942008 0.5228853
## 13 Rufous-tailed Hummingbird   dcam 0.2671073 0.10285164 0.5073230
## 14      Speckled Hummingbird   dcam 0.1924634 0.11030612 0.3612454
## 15    Stripe-throated Hermit   dcam 0.2292337 0.14697289 0.3748659
## 16      Tawny-bellied Hermit   dcam 0.2967347 0.18834135 0.4335916
## 17       Violet-tailed Sylph   dcam 0.2612181 0.15820646 0.4405280
## 18  Wedge-billed Hummingbird   dcam 0.1392272 0.05540372 0.2337274
## 19    White-whiskered Hermit   dcam 0.2601769 0.17198555 0.3587414

Discrepancy

The goodness of fit is a measured as chi-squared. The expected value for each day is the detection rate * the estimate intensity of interactions. The expected value is compared to the observed value of the actual data. In addition, a replicate dataset is generated from the posterior predicted intensity. Better fitting models will have lower discrepancy values and be Better fitting models are smaller values and closer to the 1:1 line. A perfect model would be 0 discrepancy. This is unrealsitic given the stochasticity in the sampling processes. Rather, its better to focus on relative discrepancy. In addition, a model with 0 discrepancy would likely be seriously overfit and have little to no predictive power.

fitstat<-pars_detect[pars_detect$par %in% c("fit","fitnew"),]
fitstat<-dcast(fitstat,Draw+Chain~par,value.var="estimate")

ymin<-round(min(fitstat$fit))
ymax<-round(max(fitstat$fit))
ab<-data.frame(x=0:ymax,y=0:ymax)
disc_obs<-ggplot(fitstat,aes(x=fit,y=fitnew)) + geom_point() + theme_bw() + labs(x="Discrepancy of observed data",y="Discrepancy of replicated data",col="Model")  + ggtitle("Posterior predictive check") + geom_line(data=ab,aes(x=x,y=y)) + coord_fixed() + ylim(ymin=0,ymax=max(max(c(fitstat$fit,fitstat$fitnew)))) + xlim(xmin=0,xmax=max(max(c(fitstat$fit,fitstat$fitnew))))
disc_obs

#Bayesian p-value
sum(fitstat$fitnew>fitstat$fit)/nrow(fitstat)
## [1] 0
ggsave("Figures/ObservedDiscrepancy.jpeg",width = 5,height=10)

Predicted Relationship

#Expand out pars
castdf<-dcast(pars_detect[pars_detect$par %in% c("beta1_mu","beta2_mu","beta3_mu","alpha_mu"),], Chain + Draw~par,value.var="estimate")

Posterior prediction

#Trajectories from posterior
predy<-trajF(alpha=castdf$alpha_mu,beta1=castdf$beta1_mu,trait=indat$Traitmatch,resources=indat$scaledR,beta2=castdf$beta2_mu,beta3=castdf$beta3_mu)

ggplot(data=predy,aes(x=trait)) + geom_ribbon(aes(ymin=lower,ymax=upper),alpha=0.4,fill="red")  +  theme_bw() + ylab("Interactions") + xlab("Difference between Bill and Corolla Length") + geom_point(data=indat,aes(x=Traitmatch,y=Camera)) + geom_line(aes(y=mean)) + geom_point(data=indat,aes(x=Traitmatch,y=Transect)) 

At High and Low Resource Availability

#Trajectories from posterior
predH<-trajF(alpha=castdf$alpha_mu,beta1=castdf$beta1_mu,trait=indat[indat$R=="High","Traitmatch"],resources=indat[indat$R=="High","scaledR"],beta2=castdf$beta2_mu,beta3=castdf$beta3_mu)

predM<-trajF(alpha=castdf$alpha_mu,beta1=castdf$beta1_mu,trait=indat[indat$R=="Medium","Traitmatch"],resources=indat[indat$R=="Medium","scaledR"],beta2=castdf$beta2_mu,beta3=castdf$beta3_mu)

predL<-trajF(alpha=castdf$alpha_mu,beta1=castdf$beta1_mu,trait=indat[indat$R=="Low","Traitmatch"],resources=indat[indat$R=="Low","scaledR"],beta2=castdf$beta2_mu,beta3=castdf$beta3_mu)

predhl<-melt(list(High=predH,Medium=predM,Low=predL),id.vars=colnames(predH))

colnames(predhl)[5]<-"BFlowerL"

predhl$BFlowerL<-factor(as.character(predhl$BFlowerL))
levels(predhl$BFlowerL)<-c("Low","Medium","High")

ggplot(data=predhl,aes(x=trait)) + geom_ribbon(aes(ymin=lower,ymax=upper,fill=BFlowerL),alpha=0.5)  + geom_line(aes(y=mean,col=BFlowerL),size=.8) + theme_bw() + ylab("Interactions") + xlab("Difference between Bill and Corolla Length") + geom_point(data=mindat,aes(x=Traitmatch,y=value))+ labs(fill="Resource Availability",col="Resource Availability") 

ggsave("Figures/AllRegression.jpeg",height=5,width=7)

Species Predictions

castdf<-dcast(pars_detect[pars_detect$par %in% c("beta1","beta2","beta3","alpha"),], species +Chain +Draw ~par ,value.var="estimate")

#Turn to 
castdf$species<-factor(castdf$species,levels=1:max(as.numeric(castdf$species)))

species.split<-split(castdf,list(castdf$species),drop = T)

species.traj<-list()

for(d in 1:length(species.split)){
  x<-species.split[[d]]
  index<-unique(x$species)
  
  #get data for those species
  billd<-indat[indat$jBird %in% index,]

  #scale resources
  species.traj[[d]]<-trajF(alpha=x$alpha,beta1=x$beta1,beta2=x$beta2,beta3=x$beta3,resources=billd$scaledR,trait=billd$Traitmatch)
  }

names(species.traj)<-names(species.split)

species.traj<-melt(species.traj,id.var=colnames(species.traj[[1]]))

#split out names and model
species.traj[,c("Index")]<-colsplit(species.traj$L1,"\\.",c("Index"))

spe<-merge(species.traj,jagsIndexBird,by.x="Index",by.y="jBird")

#plot and compare to original data
ggplot(data=spe,aes(x=trait)) + geom_ribbon(aes(ymin=lower,ymax=upper),alpha=0.2,fill='red') + geom_line(aes(y=mean),size=.5) + theme_bw() + ylab("Daily Interaction Rate")+ xlab("Difference between Bill and Corolla Length") + facet_wrap(~Hummingbird,scales="free",ncol=4) + geom_point(data=mindat,aes(x=Traitmatch,y=value),size=2.5) 

Species Predictions: High and Low Availability

castdf<-dcast(pars_detect[pars_detect$par %in% c("beta1","beta2","beta3","alpha"),], species +Chain +Draw ~par ,value.var="estimate")

#Turn to 
castdf$species<-factor(castdf$species,levels=1:max(as.numeric(castdf$species)))

species.split<-split(castdf,list(castdf$species),drop = T)

species.traj<-list()

for(d in 1:length(species.split)){
  x<-species.split[[d]]
  index<-unique(x$species)
  
  #get data for those species
  billd<-indat[indat$jBird %in% index,]

  sh<-trajF(alpha=x$alpha,beta1=x$beta1,beta2=x$beta2,beta3=x$beta3,resources=billd[billd$R=="High","scaledR"],trait=billd[billd$R=="High","Traitmatch"])
  
  sm<-trajF(alpha=x$alpha,beta1=x$beta1,beta2=x$beta2,beta3=x$beta3,resources=billd[billd$R=="Medium","scaledR"],trait=billd[billd$R=="Medium","Traitmatch"])

  sl<-trajF(alpha=x$alpha,beta1=x$beta1,beta2=x$beta2,beta3=x$beta3,resources=billd[billd$R=="Low","scaledR"],trait=billd[billd$R=="Low","Traitmatch"])
  
  sm<-melt(list(High=sh,Medium=sm,Low=sl),id.vars=colnames(sl))
  colnames(sm)[5]<-"Resources"
  species.traj[[d]]<-sm
  }

names(species.traj)<-names(species.split)

species.traj<-melt(species.traj,id.var=colnames(species.traj[[1]]))

#split out names and model
species.traj[,c("Index")]<-colsplit(species.traj$L1,"\\.",c("Index"))

spe<-merge(species.traj,jagsIndexBird,by.x="Index",by.y="jBird")

#plot and compare to original data
indat$rplot<-as.factor(indat$scaledR)
levels(indat$rplot)<-c("Low","Medium","High")
ggplot(data=spe,aes(x=trait)) + geom_ribbon(aes(ymin=lower,ymax=upper,fill=Resources),alpha=0.2) + geom_line(aes(y=mean,col=Resources),size=.5) + theme_bw() + ylab("Daily Interaction Rate")+ xlab("Difference between Bill and Corolla Length") + facet_wrap(~Hummingbird,scales="free",ncol=4) + geom_point(data=indat,aes(x=Traitmatch,y=Yobs,col=rplot),size=4,alpha=.75)

ggsave("Figures/SpeciesRegression.jpeg",height=6,width=7)

Species Level Interaction

castdf<-dcast(pars_detect[pars_detect$par %in% c("beta1","beta2","beta3","alpha"),], species +Chain + Draw~par,value.var="estimate")

#Turn to 
castdf$species<-factor(castdf$species,levels=1:max(as.numeric(castdf$species)))

species.split<-split(castdf,list(castdf$species),drop = T)

species.traj<-list()

for(d in 1:length(species.split)){
  dat<-species.split[[d]]
  index<-unique(dat$species)
  
  #get data for those species
  billd<-mindat[mindat$jBird %in% index,]

  #Calculate interaction effect
  species.traj[[d]]<-intF(alpha=dat$alpha,beta1=dat$beta1,x=billd[billd$value > 0 & !is.na(billd$value),'Traitmatch'],resources=billd[billd$value > 0 & !is.na(billd$value),'scaledR'],beta2=dat$beta2,beta3=dat$beta3)
  }

names(species.traj)<-names(species.split)
species.traj<-melt(species.traj,id.var=colnames(species.traj[[1]]))

#split out names and model
species.traj[,c("Index")]<-colsplit(species.traj$L1,"\\.",c("Index"))

spe<-merge(species.traj,jagsIndexBird,by.x="Index",by.y="jBird")

#match colnames

#plot and compare to original data
ggplot(data=spe,aes(x=x)) + geom_ribbon(aes(ymin=lower,ymax=upper,fill=Hummingbird),alpha=0.3) + geom_line(aes(y=mean,col=Hummingbird),size=1) + theme_bw() + xlab("Difference between Bill and Corolla Length")  + ylab("Effect of Resources on Trait Difference") + facet_wrap(~Hummingbird,scales="free",ncol=3)
ggsave("Figures/SpeciesInteraction.jpeg",height=6,width=7)

These plots can be tricky to interpret if one forgets that trait matching as a covariate is a distance. Therefore remember that a positive slope in the plot above indiciates, “As resources increase species use flowers less similiar to their bill lengths”.

Which species did we predict well?

smat<-pars_detect %>% filter(par=="E") %>% group_by(species) %>% summarize(fit=sum(estimate)) %>% arrange(desc(fit))
smat<-merge(smat,jagsIndexBird,by.x="species",by.y="jBird")
smat
##    species        fit               Hummingbird
## 1        1   4198.651            Andean Emerald
## 2        2  90919.540        Booted Racket-tail
## 3        3 209789.915                Brown Inca
## 4        4  38912.967       Buff-tailed Coronet
## 5        5  92914.052             Collared Inca
## 6        6  24295.666         Crowned Woodnymph
## 7        7  18216.200   Fawn-breasted Brilliant
## 8        8  31442.695         Gorgeted Sunangel
## 9        9  13621.877   Green-crowned Brilliant
## 10      10  41670.690   Green-fronted Lancebill
## 11      11  12204.590             Hoary Puffleg
## 12      12  12602.366    Purple-bibbed Whitetip
## 13      13    847.226 Rufous-tailed Hummingbird
## 14      14  95813.982      Speckled Hummingbird
## 15      15  94792.828    Stripe-throated Hermit
## 16      16 151151.034      Tawny-bellied Hermit
## 17      17 281231.597       Violet-tailed Sylph
## 18      18  33752.496  Wedge-billed Hummingbird
## 19      19  83972.069    White-whiskered Hermit

Discrepancy matrix

dmat<-pars_detect %>% filter(par=="E") %>% group_by(species,plant) %>% summarize(fit=sum(estimate)) %>% arrange(desc(fit))
dmat<-merge(dmat,jagsIndexBird,by.x="species",by.y="jBird")
dmat<-merge(dmat,jagsIndexPlants,by.x="plant",by.y="jPlant")
head(dmat,10)
##    plant species        fit              Hummingbird
## 1      1       3 6192.83019               Brown Inca
## 2      1       2  124.05703       Booted Racket-tail
## 3      1      18   36.06558 Wedge-billed Hummingbird
## 4      1       9   50.77790  Green-crowned Brilliant
## 5      1       4   29.23103      Buff-tailed Coronet
## 6      1      10  259.95336  Green-fronted Lancebill
## 7      1      17 2121.24617      Violet-tailed Sylph
## 8      1       7  822.22259  Fawn-breasted Brilliant
## 9      1      16 1712.61751     Tawny-bellied Hermit
## 10     1       5  536.80132            Collared Inca
##                Iplant_Double
## 1  Alloplectus tetragonoides
## 2  Alloplectus tetragonoides
## 3  Alloplectus tetragonoides
## 4  Alloplectus tetragonoides
## 5  Alloplectus tetragonoides
## 6  Alloplectus tetragonoides
## 7  Alloplectus tetragonoides
## 8  Alloplectus tetragonoides
## 9  Alloplectus tetragonoides
## 10 Alloplectus tetragonoides
#get order
hord<-dmat %>% group_by(Hummingbird) %>% summarize(n=sum(fit)) %>% arrange(desc(n)) %>% .$Hummingbird
pord<-dmat %>% group_by(Iplant_Double) %>% summarize(n=sum(fit)) %>% arrange(n) %>%  .$Iplant_Double
dmat$Hummingbird<-factor(dmat$Hummingbird,levels=hord)
dmat$Iplant_Double<-factor(dmat$Iplant_Double,levels=pord)
ggplot(dmat,aes(x=Hummingbird,Iplant_Double,fill=fit)) + geom_tile() + theme_bw() + scale_fill_continuous("Discrepancy",low="blue",high="red")

Resource Abundance Functions

Let’s take a closer look at distribution of interaction effect posteriors values for each species.

post<-pars_detect %>% filter(par %in% "beta2") %>% group_by(Hummingbird) %>% summarize(mean=mean(estimate),median=median(estimate),lower=quantile(probs=0.025,estimate),upper=quantile(probs=0.975,estimate)) %>% melt(id.vars='Hummingbird')
ggplot(pars_detect[pars_detect$par %in% "beta2",],aes(x=estimate)) + geom_histogram() + facet_wrap(~Hummingbird,scales='free',ncol=4) + geom_vline(data=post,aes(xintercept=value,col=variable))

Interaction density functions

Let’s take a closer look at distribution of interaction effect posteriors values for each species.

post<-pars_detect %>% filter(par %in% "beta3") %>% group_by(Hummingbird) %>% summarize(mean=mean(estimate),median=median(estimate),lower=quantile(probs=0.025,estimate),upper=quantile(probs=0.975,estimate)) %>% melt(id.vars='Hummingbird')
ggplot(pars_detect[pars_detect$par %in% "beta3",],aes(x=estimate)) + geom_histogram() + facet_wrap(~Hummingbird,scales='free',ncol=4) + geom_vline(data=post,aes(xintercept=value,col=variable))

Trait-matching and Bill Length

Do species with long bill lengths have positive traitmatching effects?

#species names
b<-pars_detect[pars_detect$par %in% "beta1",]

#traits
b<-merge(b,hum.morph,by.x="Hummingbird",by.y="English")

post<-b %>% filter(par %in% "beta1") %>% group_by(Hummingbird) %>% summarize(mean=mean(estimate),median=median(estimate),lower=quantile(probs=0.025,estimate),upper=quantile(probs=0.975,estimate),quantile_l=quantile(estimate)[[1]],quantile_u=quantile(estimate)[[2]]) %>% melt(id.vars='Hummingbird')

#get order of mean posterior
ord<-post %>% filter(variable=="mean") %>% arrange(value) %>% .$Hummingbird

b$Hummingbird<-factor(b$Hummingbird,levels=ord)
ggplot(b,aes(y=estimate,x=as.factor(Total_Culmen),group=Hummingbird)) + geom_violin(fill='grey50') + coord_flip()  + ggtitle("Trait-matching and Bill Length") + theme_bw()

Interaction and Bill Length

Do species with long bill lengths have positive interaction effects?

#species names
b<-pars_detect[pars_detect$par %in% "beta3",]

#traits
b<-merge(b,hum.morph,by.x="Hummingbird",by.y="English")

post<-b %>% filter(par %in% "beta3") %>% group_by(Hummingbird) %>% summarize(mean=mean(estimate),median=median(estimate),lower=quantile(probs=0.025,estimate),upper=quantile(probs=0.975,estimate),quantile_l=quantile(estimate)[[1]],quantile_u=quantile(estimate)[[2]]) %>% melt(id.vars='Hummingbird')

#get order of mean posterior
ord<-post %>% filter(variable=="mean") %>% arrange(value) %>% .$Hummingbird

b$Hummingbird<-factor(b$Hummingbird,levels=ord)
ggplot(b,aes(y=estimate,x=Hummingbird,fill=Total_Culmen)) + geom_violin() + coord_flip() + scale_fill_continuous(low='blue',high='red') + ggtitle("Interaction Effect and Bill Length") + theme_bw()

Estimated niche breadth

castdf<-dcast(pars_detect[pars_detect$par %in% c("beta1","beta2","beta3","alpha"),], species +Chain + Draw~par,value.var="estimate")

#Turn to 
castdf$species<-factor(castdf$species,levels=1:max(as.numeric(castdf$species)))

species.split<-split(castdf,list(castdf$species),drop = T)

species.traj<-lapply(species.split,function(dat){
  index<-unique(dat$species)
  
  #get data for those species
  billd<-indat[indat$jBird %in% index,]
  
  d<-data.frame(alpha=dat$alpha,beta1=dat$beta1,beta2=dat$beta2,beta3=dat$beta3)
  
  #fit regression for each input estimate
  sampletraj<-list()
  
  for (y in 1:nrow(d)){
    v=exp(d$alpha[y] + d$beta1[y] * billd$Traitmatch + d$beta2[y] * billd$scaledR + d$beta3[y] * billd$Traitmatch*billd$scaledR)
    
    sampletraj[[y]]<-data.frame(x=as.numeric(billd$Traitmatch),y=as.numeric(v),r=as.numeric(billd$scaledR),jBird=billd$jBird,jPlant=billd$jPlant,jTime=billd$jTime)
  }
  
  sample_all<-rbind_all(sampletraj)
})
  
species.traj<-rbind_all(species.traj)

Mean Estimates for Corolla Sizes

species.mean<-species.traj %>% group_by(jBird,jPlant,r) %>% summarize(Traitmatch=unique(x),phi=mean(y))

tomerge<-indat %>% select(jBird,jPlant,Hummingbird,Iplant_Double) %>% distinct()

species.mean<-merge(species.mean,tomerge)

#get corolla sizes
species.mean<-merge(species.mean,fl.morph,by.x="Iplant_Double", by.y="Group.1")

#bill order
ord<-hum.morph %>% arrange(Total_Culmen) %>% .$English
species.mean$Hummingbird<-factor(species.mean$Hummingbird,levels=ord)

#add level to hum.morph to match naming convention
species.mean<-merge(species.mean,hum.morph[,c("English","Total_Culmen")],by.x="Hummingbird",by.y="English")

p<-ggplot(species.mean,aes(x=Corolla,y=phi,col=as.factor(r))) + geom_line(size=.9) + geom_vline(aes(xintercept=Total_Culmen),linetype='dashed') + facet_wrap(~Hummingbird,ncol=3,scales="free_y")  + theme_bw() + ylab("Probability of Interaction") + scale_color_manual("Resource Availability",labels=c("Low","Medium","High"),values=c("Black","Blue","Red")) + xlab("Flower Corolla Length (mm)") 
p

ggsave("Figures/ResponseCurves.jpeg",height=6.5,width=9)

Niche Breadth

species.mean<-species.traj %>% group_by(jBird,jPlant,r) %>% summarize(Traitmatch=unique(x),phi=mean(y),phi_low=quantile(y,0.05),phi_high=quantile(y,0.95))

#merge names
species.mean<-merge(species.mean,jagsIndexBird)
species.mean<-merge(species.mean,jagsIndexPlants)

#get corolla sizes
species.mean<-merge(species.mean,fl.morph,by.x="Iplant_Double", by.y="Group.1")

#bill order
ord<-hum.morph %>% arrange(Total_Culmen) %>% .$English
species.mean$Hummingbird<-factor(species.mean$Hummingbird,levels=ord)

#add level to hum.morph to match naming convention
species.mean<-merge(species.mean,hum.morph[,c("English","Total_Culmen")],by.x="Hummingbird",by.y="English")

#label factor
species.mean$r<-as.factor(species.mean$r)
levels(species.mean$r)<-c("Low","Medium","High")
ggplot(species.mean) + geom_ribbon(alpha=0.7,aes(x=Corolla,ymin=phi_low,ymax=phi_high,fill=r)) + theme_bw() + facet_wrap(~Hummingbird,scales="free",ncol=4)+ ggtitle("Niche Breadth") + geom_vline(aes(xintercept=Total_Culmen),linetype='dashed') + geom_line(aes(x=Corolla,y=phi,fill=r)) + ylab("Probability of Interaction") + xlab("Corolla Length (mm)") + scale_fill_manual("Resource Availability",values=c("Grey70","Grey20","Black"))

ggsave("Figures/NicheBreadth.jpeg",height=8,width=9)

Range plots

nsplit<-split(species.mean,species.mean$r)
makeR<-function(x){
  
  #input matrix
  aggm<-matrix(nrow=nrow(jagsIndexBird),ncol=nrow(jagsIndexPlants),data=0)
  for (j in 1:nrow(x)){
    aggm[x[j,"jBird"],x[j,"jPlant"]]<-rpois(1,lambda=x[j,"phi"])
  }
  
  aggm<-melt(aggm)
  colnames(aggm)<-c("jBird","jPlant","P")
  tomerge<-species.mean %>% select(jBird,jPlant,Corolla,Hummingbird,Traitmatch) %>% distinct()
  aggm<-merge(aggm,tomerge)
  return(aggm)
}

nstat<-lapply(nsplit,function(x){
  netstat<-melt(lapply(1:50,function(k) makeR(x)),id.vars=c("jBird","jPlant","Hummingbird","Traitmatch","Corolla","P"))
  colnames(netstat)[7]<-"Iteration"
  return(netstat)
})

names(nstat)<-c("Low","Medium","High")
nstat<-melt(nstat,colnames(nstat[[1]]))

Predicted density

ggplot(nstat[nstat$P==1,],aes(x=Corolla,fill=L1)) + geom_density(alpha=0.6) + facet_wrap(~Hummingbird,scales='free',nrow=5) + scale_fill_discrete("Resource Availability") + theme_bw()

rangestat<-nstat %>% filter(P==1) %>% group_by(Hummingbird,L1) %>% summarize(mean=mean(Corolla),lower=quantile(Corolla,0.05),upper=quantile(Corolla,0.95))
  
ggplot(rangestat,aes(x=Hummingbird,ymin=lower,ymax=upper,ymean=mean,col=L1)) + geom_linerange(alpha=0.5,position=position_dodge(width=.75),size=2) + geom_point(aes(y=mean),position=position_dodge(width=.75)) + labs(col="Resource Availability",y="Corolla Length (mm)") + theme_bw() + coord_flip() + theme(legend.position="bottom")

Predicted Mean

meanstat<-nstat %>% filter(P==1) %>% group_by(Hummingbird,L1,Iteration) %>% summarize(a=mean(Corolla))%>% summarize(mean=mean(a),lower=quantile(a,0.05),upper=quantile(a,0.95))
  
ggplot(meanstat,aes(x=Hummingbird,ymin=lower,ymax=upper,ymean=mean,col=L1)) + geom_linerange(alpha=0.5,position=position_dodge(width=.6),size=2) + geom_point(aes(y=mean),position=position_dodge(width=.6)) + labs(col="Resource Availability",y="Corolla Length (mm)") + theme_bw() + coord_flip() + theme(legend.position="bottom")

Predicted Min

meanstat<-nstat %>% filter(P==1) %>% group_by(Hummingbird,L1,Iteration) %>% summarize(a=min(Corolla))%>% summarize(mean=mean(a),lower=quantile(a,0.05),upper=quantile(a,0.95))
  
ggplot(meanstat,aes(x=Hummingbird,ymin=lower,ymax=upper,ymean=mean,col=L1)) + geom_linerange(alpha=0.5,position=position_dodge(width=.75),size=2) + geom_point(aes(y=mean),position=position_dodge(width=.75)) + labs(col="Resource Availability",y="Corolla Length (mm)") + theme_bw() + coord_flip() + theme(legend.position="bottom")

Predicted Max

meanstat<-nstat %>% filter(P==1) %>% group_by(Hummingbird,L1,Iteration) %>% summarize(a=min(Corolla))%>% summarize(mean=mean(a),lower=quantile(a,0.05),upper=quantile(a,0.95))
  
ggplot(meanstat,aes(x=Hummingbird,ymin=lower,ymax=upper,ymean=mean,col=L1)) + geom_linerange(alpha=0.5,position=position_dodge(width=.6),size=2) + geom_point(aes(y=mean),position=position_dodge(width=.6)) + labs(col="Resource Availability",y="Corolla Length (mm)") + theme_bw() + coord_flip() + theme(legend.position="bottom")

Generate network

#Split by resource
nsplit<-split(species.mean,species.mean$r)

makeN<-function(x){
  
  #input matrix
  aggm<-matrix(nrow=nrow(jagsIndexBird),ncol=nrow(jagsIndexPlants),data=0)
  for (j in 1:nrow(x)){
    aggm[x[j,"jBird"],x[j,"jPlant"]]<-rpois(1,lambda=x[j,"phi"])
  }
  #calculate network statistic
  nstat<-networklevel(aggm,index=c("weighted connectance","weighted NODF","niche overlap"),level="lower")
}

nstat<-lapply(nsplit,function(x){
  netstat<-melt(t(sapply(1:500,function(k) makeN(x)))) 
  colnames(netstat)<-c("Iteration","Metric","value")
  return(netstat)
})

names(nstat)<-c("Low","Medium","High")
nstat<-melt(nstat,colnames(nstat[[1]]))

ggplot(nstat,aes(x=value,fill=L1)) + geom_density(alpha=0.6) + facet_wrap(~Metric,scales='free',nrow=3) + scale_fill_discrete("Resource Availability") + theme_bw()

ggsave("Figures/NetworkStatistics.jpeg",height=5,width=6,dpi=600)

Compared to raw visits

sindat<-split(indat,indat$scaledR)
ndat<-lapply(sindat,function(x){
  web<-acast(x[x$Yobs==1,],Hummingbird~Iplant_Double,value.var = "Yobs", fun= function(x){ length(unique(x))})
  nstat<-t(networklevel(web,index=c("weighted connectance","weighted NODF","niche overlap"),level="lower"))
})
names(ndat)<-c("Low","Medium","High")
ndat<-melt(ndat)

colnames(ndat)<-c("Var1","Metric","value","L1")

ggplot(nstat,aes(x=value,fill=L1)) + geom_density(alpha=0.6) + facet_wrap(~Metric,scales='free',nrow=2) + scale_fill_discrete("Resource Availability") + geom_vline(data=ndat,aes(xintercept=value,col=L1),linetype="dashed") + scale_color_discrete(guide='none')

save.image("ObservedModel.RData")